
options(error=stop)
rm(list=ls())

load("./RDataFiles/results1.RData")

Results = result

### Create a summary table 

measure = c("Mean","Bias","SE","RMSE","Cov","OutofRange")
results.summary = array(NA,c(dim(Results)[2:7],length(measure)))
dimnames.results.summary = dimnames(Results)[2:7]
dimnames.results.summary[[7]] = measure
dimnames(results.summary) = dimnames.results.summary

sd      = function(x) sqrt(var(x,na.rm=TRUE))
rmse    = function(x) sqrt(mean((x-mean(x))^2))
outofrange = function(x) sum(abs(x)>1)
cov     = function(x) {
    se = sd(x)
    return(mean( abs( (x-Delta.true)/se ) < qnorm(0.975,0,1) ) )
}

results.summary[,,,,,,"Mean"] = apply(Results,2:7,mean,na.rm=TRUE)
results.summary[,,,,,,"Bias"] = results.summary[,,,,,,"Mean"] - Delta.true
results.summary[,,,,,,"SE"]   = apply(Results,2:7,sd)/sqrt(rep.times)
results.summary[,,,,,,"Cov"]   = apply(Results,2:7,cov)
results.summary[,,,,,,"RMSE"]       = apply(Results,2:7,rmse)
results.summary[,,,,,,"OutofRange"] = apply(Results,2:7,outofrange)


### Creating Table 2 in the main text

scenarios = c("All correct", "M1 correct", "M2 correct", "M3 correct", "All wrong")
report = array(NA, dim=c(5,dim(results.summary)[5:7]),
               dimnames = c(list(scenario = scenarios), 
                            dimnames(results.summary)[5:7]))
report[1,,,] = results.summary[1,1,1,1,,,]
report[2,,,] = results.summary[2,1,1,1,,,]
report[3,,,] = results.summary[1,1,2,2,,,]
report[4,,,] = results.summary[1,2,1,2,,,]
report[5,,,] = results.summary[2,2,2,2,,,]

library(xtable)
n.index = "n=500"
est.index = c("reg","b-ipw","g","tr","b-tr")
table.mean = format(round(report[,n.index,est.index,"Bias"],3),3)
table.se   = format(round(report[,n.index,est.index,"SE"],3),3)
table.mean.se = matrix(paste(table.mean, "(", table.se, ")",sep=""),
                       dim(table.mean))
dimnames(table.mean.se) = dimnames(table.mean)
table.mean.se

xtable(table.mean.se)
xtable(format(round(report[,n.index,est.index,"RMSE"],3),3))



### Creating Figure 2 in the main text

oor.tr = apply(Results[,1,2,1,2,"n=500",],1,function(x){
    if(abs(x["tr"])>1) return(x["b-tr"])
    return(NA)
})
oor.tr = c(oor.tr);   oor.tr = oor.tr[!is.na(oor.tr)]
hist(oor.tr,20,prob=TRUE)
lines(density(oor.tr, adjust=2))
abline(v=Delta.true,col="red",lwd=2)

temp = Results[,1,2,1,2,"n=500","tr"]
hist(temp[abs(temp)<50])


library(ggplot2)
pdf("./Figures/b-tr-estimates.pdf",onefile=FALSE, width = 6.0, height = 4.3)
plot.data = data.frame(oor.tr = oor.tr, true = Delta.true)
ggplot(data = plot.data, aes(oor.tr)) + 
    theme_bw(base_size = 12, base_family = "") +
    geom_histogram(binwidth = 0.015, aes(y = ..density..),
                   colour = "darkgreen", fill = "white") + 
    geom_density(adjust = 2, lwd = 0.75) + 
    geom_vline(xintercept = Delta.true, colour = "red", lwd = 1) +
    labs(x="b-tr estimates") + labs(y="density") +
    theme(axis.title.y = element_text(size = rel(1.5)),
          axis.title.x = element_text(size = rel(1.5)),
          axis.text = element_text(size=rel(1.2)))
dev.off()
    

### Creating boxplots in Table 4 in the Appendix

source("liteep.R"); source("1.6.1 Boxplot.R")
paras = list(
    pi_true = c(1,2),       # 1 = correct, 3 = incorrect
    deltad_true = c(1,2),   # 1 = correct, 2 = incorrect
    delta_true = c(1,2),    # 1 = correct, 2 = incorrect
    op_true = c(1,2),       # 1 = correct, 2 = incorrect
    n.cand = c(500)
)
para.values = ParaValues(paras)

apply(para.values,1:5,function(x) Boxplot(x[5],x[1:4]))

Boxplot(500,c(1,1,1,1))
















